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ABSTRACT 


Two  possible  mechanisms  for  the  accidental  ignition  of  liquid  propellants  are  studied:  (i)  The  heating  of  the 
liquid  surrounding  a  strongly  compressed  gas  bubble,  and  (it)  The  heating  due  to  viscous  dissipation  in  a 
rapidly  compressed  drop.  Contrary  to  previous  work,  for  the  first  problem  the  heating  of  the  gas  is  calculated 
precisely  rather  than  estimated  on  the  basis  of  an  adiabatic  model,  while  for  the  second  one  proper  allowance 

is  made  for  the  temperature  dependence  of  the  liquid  viscosity.  ,  , 

For  the  first  problem  the  compression  due  to  a  pressure  pulse  and  a  pressure  step  have  been  studied. 
In  the  first  case,  the  maximum  response  is  found  when  the  pulse  time  scale  is  - comparabk  to  ' that  for  the 
collapse  of  the  bubble.  The  temperature  reached  by  the  liquid  can  be  m  excess  of  1,000  K,  which  presuma  y 
would  lead  to  vaporization  not  accounted  for  in  the  model.  The  duration  of  such  high  temperatures  however 
is  probably  too  limited  for  direct  ignition.  However,  if  the  propellant  vaporizes,  its  vapor  can  conceivably 
be  heated  during  the  next  compression  cycle  of  the  bubble  oscillatory  motion.  Some  limited  results  for  the 

case  of  non-spherical  collapse  of  the  bubble  are  also  presented. 

For  the  viscous  heating  problem  a  marked  effect  of  the  temperature  dependence  of  viscosity  has  been 
found.  Furthermore,  although  very  high  temperatures  develop  when  the  velocity  with  which  the  <kop  is 
compressed  is  held  constant,  this  only  happens  when  the  pressure  is  so  large  that  the  assumption  of  a 
constant  velocity  is  unrealistic.  The  effect  of  a  wall  slip  boundary  condition  has  also  been  studied  and  fo 
to  be  large.  The  results  indicate  that  the  process  under  consideration  is  strongly  influenced  by  subtle  aspects 

that  cannot  be  ignored  even  as  a  first  approximation.  ,  .  „  m 

In  general,  all  the  results  indicate  that  ignition  is  in  principle  possible,  although  marginally  so.  While 
this  finding  is  generally  consistent  with  experiment,  it  may  also  explain  why  it  is  difficult  to  deduce  definite 
conclusions  from  the  available  experimental  evidence  with  its  apparently  inconsistent  resu  ts. 
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Id?hOM  AND  SHOULD  ■ NOT  BE  CONSTRUED  AS  AN  OFFICIAL  DEPARTMENT  OF  THE  ARMY 
PWlS  POLICY  OR  DEC, SION,  UNLESS  SO  DESIGNATED  BY  OTHER  DOCUMENTATION 
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Fig.  1.  Normalized  radius  R/Ro  versus  the  normalized  time  t/tR  [where  Rq  =  1  mm  is  the  initial  radius 
and  tR  is  the  Rayleigh  collapse  time  defined  in  (2.3.2)]  for  an  air  bubble  subject  to  Gaussian  compression 
pulse  of  width  tp.  The  dashed  line  is  for  tR/tp  =  0.1,  the  solid  line  for  iR/tp  =  1,  and  the  dotted  line  for 
tR/tp  =  10.  The  initial  pressure  is  1  bar,  the  maximum  overpressure  10  bars,  and  tR  ~  29  ps. 

Fig.  2.  Temperature  attained  by  the  gas  at  the  bubble  center  (squares)  and  by  the  liquid  at  the  bubble 
surface  (circles)  at  the  end  of  the  first  collapse  for  tR/tp  =  0.1,  0.2,  0.5,  1,  2,  5,  and  10.  The  dotted  line  is 
the  adiabatic  prediction.  The  collapse  is  caused  by  a  Gaussian  pressure  pulse  with  width  tp  and  maximum 
amplitude  10  bars.  The  initial  bubble  radius  is  1  mm  and  the  initial  ambient  pressure  1  bar. 

Fig.  3.  Center  gas  temperature  for  the  three  cases  of  Fig.  1  as  a  function  of  t/tR.  The  dashed  curve  is  for 
tR/tp  =  0.1,  the  solid  one  for  tR/tp  =  1,  and  the  dotted  one  for  tR/tp  =  10. 

pjg,  4.  Liquid  surface  temperature  for  the  three  cases  of  Fig.  1  as  a  function  of  t/tR.  The  dashed  curve  is 
for  tR/tp  =  0.1,  the  solid  one  for  tR/tp  =  1,  and  the  dotted  one  for  tR/tP  =  10. 

Fig.  5.  Wall  heat  fluxes  for  the  three  cases  of  Fig.  1  as  a  function  of  t/tR.  The  dashed  curve  is  for  tR/tp 
-  0.1,  the  solid  one  for  tR/tP  =  1,  and  the  dotted  one  for  tR/tp  =  10. 

Fig.  6.  Amount  of  heat  delivered  to  the  liquid  (2.3.4)  as  a  function  of  time  for  the  three  cases  of  Fig.  1. 
The  dashed  curve  is  for  tR/tp  =  0.1,  the  solid  one  for  tR/tp  =  1,  and  the  dotted  one  for  tR/tp  -  10. 

Fig.  7.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  tR  -  tp,  Ap  =  10  bars,  and  R0  = 
0.1  mm.  The  Rayleigh  time  is  tR  =  2.9  ps. 

Fig.  8.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  tR  —  tp,  A p  =  10  bars,  and  Ro  =  10 
mm.  The  Rayleigh  time  is  tR  =  0.29  ms. 

Fig.  9.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  tR  =tp,  Ap  =  100  bars,  and  Ro  = 
0.1  mm.  The  Rayleigh  time  is  tR  -  0.92  ps. 

Fig.  10.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  tR  =  tp,  Ap  =  100  bars,  and  R0  = 

1  mm.  The  Rayleigh  time  is  tR  =  9.2  ps. 

Fig.  11.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  tR  —  tpi  Ap  =  100  bars,  and  Ro  = 
10  mm.  The  Rayleigh  time  is  tR  =  92  ps. 

Fig.  12.  Normalized  radius  R/Rq  versus  the  normalized  time  t/tR  for  an  air  bubble  with  initial  radius  Ro 
=  1  mm  subject  to  a  step  pressure  increase  Ap  =  10  bars.  The  Rayleigh  collapse  time  is  tR  =  29  ps. 

Fig.  13.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  Ap  = 
10  bars  and  R0  =  0.1  mm.  The  Rayleigh  time  is  tR  =  2.9  ps. 

Fig.  14.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  Ap  = 
10  bars  and  /20  =  1  mm.  The  Rayleigh  time  is  tR  =  29  ps. 

Fig.  15.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  Ap  = 
10  bars  and  R0  =  10  mm.  The  Rayleigh  time  is  tR  =  0.29  ms. 
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Fig.  16.  Normalized  radius  R/Ro  versus  the  normalized  time  t/tR  for  an  air  bubble  with  initial  radius  Rq 
—  1  mm  subject  to  a  step  pressure  increase  A p  =  100  bars.  The  Rayleigh  collapse  time  is  tR  =  9.2  /is. 

Fig.  17.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  A p  = 
100  bars  and  R0  =  0.1  mm.  The  Rayleigh  time  is  tR  =  0.92  ps. 

Fig.  18.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  A p  = 
100  bars  and  R0  =  1  mm.  The  Rayleigh  time  is  tR  =  9.2  ps. 

Fig.  19.  Core  gas  (dashed)  and  bubble  surface  (solid)  temperatures  for  a  step  pressure  increase  with  A p  = 
100  bars  and  R0  =  10  mm.  The  Rayleigh  time  is  tR  =  92  ps. 

Fig.  20.  Bubble  surface  temperature  for  a  step  pressure  increase  with  A p  =  100  bars  and  R0  =  10  mm. 
The  Rayleigh  time  is  tR  =  92  ps.  These  results  differ  from  those  of  Fig.  19  because  here  the  liquid  thermal 
conductivity  is  0.075  rather  than  0.15  W/mK  and  the  gas  adiabatic  index  is  5/3  rather  than  7/5. 

Fig.  21.  Example  of  a  numerically  generated  orthogonal  grid  for  the  solution  of  the  Laplace  and  energy 
equation  in  the  gas.  Here  the  initial  bubble  has  a  radius  Rq  =  1  mm,  with  a  rigid  wall  1.2  mm  below  the 
initial  bubble  center.  The  collapse  is  triggered  by  a  pressure  step  Ap  =  10  bars.  The  grid  is  shown  at  the 
time  t  =  26.64  ps.  The  Rayleigh  time  is  tR  =  29  ps. 

Fig.  22.  The  left  figure  shows  isotherms  calculated  with  the  full  model  for  Rq  —  1  mm,  with  a  rigid  wall 
1.2  mm  below  the  initial  bubble  center  and  a  pressure  step  Ap  =  10  bars  at  t  =  26.64  ps.  The  right  figure 
shows  the  corresponding  lines  of  constant  $.  Here  the  Rayleigh  time  is  Ir  =  29  ps. 

Fig.  23.  The  gas  temperature  along  the  axis  of  symmetry  of  the  bubble  (from  the  closest  to  the  farthest 
point  from  the  wall)  for  the  case  of  the  previous  figure. 

Fig.  24.  Comparison  between  the  gas  temperatures  predicted  by  the  full  model  (dashed  line)  and  the 
boundary  layer  approximation  (solid  line)  for  the  spherical  collapse  of  an  air  bubble  with  Rq  —  1  mm  and  a 
step  pressure  increase  Ap  =  10  bars. 

Fig.  25.  Comparison  between  the  heat  fluxes  at  the  jet  tip  predicted  by  the  full  model  (dashed  line)  and 
the  boundary  layer  approximation  (solid  line)  for  an  asymmetric  collapse.  The  case  is  that  of  Figs.  21  to  23. 

Fig.  26.  Successive  bubble  configurations  for  a  collapse  induced  by  a  sudden  overpressure  with  Ap  =  10 
bars,  po  =  1  bar.  The  initial  bubble  radius  is  R0  =  1  mm  and  the  rigid  wall  is  initially  1.2  mm  below  the 
bubble  center. 

Fig.  27.  Successive  shapes  of  a  1  mm-radius  air  bubble  subjected  to  a  Gaussian  pressure  pulse  (2.3.1)  with 
with  tR/tp  =  4  and  Ap  =  10  bars.  Notice  that  the  jet  continues  to  move  toward  the  opposite  bubble  surface 
even  while  the  rest  of  the  bubble  is  rebounding. 

Fig.  28.  The  local  heat  flux  along  the  bubble  surface  for  a  step  pressure  increase  Ap  =  10  bars  as  computed 
from  the  boundary-layer  model  at  different  times  for  Ro  =  1  mm.  The  abscissa  denotes  the  node  number. 
At  each  time  step  sixty  equispaced  nodes  were  used. 

Fig.  29.  Estimated  jet  tip  liquid  temperature  as  a  function  of  time  for  the  case  of  Figs.  21  to  23  (solid 
line).  The  dashed  line  is  the  liquid  temperature  for  a  spherical  collapse. 
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Fig.  1  The  left  part  of  the  figure  is  a  sketch  of  the  experiment  that  motivates  this  study.  The  right  part 
shows  the  process  simulated  in  the  present  computations  with  the  frame  reference  centered  on  the  plane  of 
symmetry. 

Fig.  2  Distribution  of  the  radial  velocity  u  across  the  gap  at  the  edge  of  the  plates  for  a  constant  velocity  of 
approach  Vo  =  1  m/s  and  adiabatic  walls.  The  lines  are,  in  descending  order,  for  h/h0  =  0.00369,  0.00451, 
0.00551,  0.00822,  0.0123,  0.0223,  0.0497,  0.819. 

Fig.  3  Temperature  distribution  across  the  gap  at  the  edge  of  the  plates  for  a  constant  velocity  of  approach 
Vo  =  1  m/s  and  an  adiabatic  wall.  The  lines  are,  in  descending  order,  for  h/h0  =  0.00369,  0.00451,  0.00551, 
0.00822,  0.0123,  0.0223,  0.0497,  0.819. 

Fig.  4  Maximum  liquid  temperature  (occurring  on  the  plates  at  r  =  R)  as  a  function  of  the  dimensionless 
gap  width  h/h( 0)  for  the  adiabatic-wall  case  and  a  fixed  velocity  of  approach  V0  =  1  m/s.  The  solid  line  is 
for  a  temperature-dependent  viscosity,  the  dashed  line  for  a  constant  viscosity. 

Fig.  5  The  radial  velocity  on  the  plane  of  symmetry  at  the  plates’  edge  r  =  R  as  a  function  of  the  gap  width 
h/h(0)  for  a  constant  velocity  U0  =  1  m/s  and  adiabatic  walls.  The  solid  line  is  for  a  temperature-dependent 
viscosity,  the  dashed  line  for  a  constant  viscosity. 

Fig.  6  The  solid  lines  are  the  radial  velocity  on  the  plane  of  symmetry  at  the  plates’  edge  r  =  R  for 
constant  velocities  Vo  =  1  m/s  (upper  curve)  and  0.5  m/s  as  functions  of  h/h(Q).  The  dashed  lines  are  the 
corresponding  temperature  increases  at  the  plates  for  adiabatic  boundary  conditions.  Units  are  in  m/s  for 
the  velocities  and  K  for  the  temperature  increases. 

Fig.  7  Distribution  of  the  radial  velocity  across  the  gap  at  the  edge  of  the  plates  for  a  constant  velocity 
of  approach  V0  =  1  m/s  and  an  isothermal  wall.  The  lines  are  for  the  same  values  of  h/h( 0)  as  in  Fig.  2. 
Specifically,  in  descending  order,  h/ho  =  0.00369,  0.00451,  0.00551,  0.00822,  0.0123,  0.0223,  0.0497,  0.819. 

Fig.  8  Temperature  distributions  across  the  gap  at  r  —  R  for  an  isothermal  wall  for  the  same  values  of 
h/h( 0)  as  in  the  previous  figure. 

Fig.  9  The  solid  lines  are  the  radial  velocity  on  the  plane  of  symmetry  at  the  plates’  edge  r  =  R  for 
constant  velocities  Vo  =  1  m/s  (upper  curve)  and  0.5  m/s  as  functions  of  h/h{ 0).  The  dashed  lines  are  the 
corresponding  temperature  increases  at  the  plates  for  isothermal  boundary  conditions.  Units  are  in  m/s  for 
the  velocities  and  K  for  the  temperature  increases. 

Fig.  10  Same  as  the  previous  figure  for  a  variable  velocity  of  approach  controlled  by  the  inertia  of  the 
plates.  The  three  curves  are,  in  ascending  order,  for  M  =  0.2^  1,  and  2  kg  and  the  walls  are  adiabatic. 

Fig.  11  Same  as  Fig.  10  for  isothermal  walls. 

Fig.  12  Temperature  rises  in  the  presence  of  slip,  with  adiabatic  walls  and  a  fixed  velocity  Vo  1  m/s.  The 
three  curves  are,  in  ascending  order,  for  A  =  0,  10,  and  100  nm. 


1  Introduction 


The  increasing  technological  importance  of  liquid  monopropellants  such  as  LGP  1845  and  1846  motivates 
a  strong  interest  in  the  safe  operational  limits  of  their  use  (see  e.g.  Knapton  et  al.  1992  for  a  review).  In 
particular,  the  propellant  behavior  under  conditions  of  high  pressures,  strong  impact,  or  high-speed  flow  has 
a  direct  bearing  on  safety  issues.  In  view  of  its  thermal  nature,  for  ignition  to  occur,  the  mechanical  energy 
associated  with  these  processes  must  somehow  lead  to  the  formation  of  so-called  ‘hot  spots’.  A  number  of 
mechanisms  for  this  to  happen  have  been  described  in  the  literature  (see  e.g.  Field  et  al.  1982;  Field  et  al 
1992).  In  the  work  described  in  this  report,  following  our  original  proposal,  we  focus  on  two: 

1.  The  heating  of  the  liquid  surrounding  a  strongly  compressed  gas  bubble  avoiding  the  assumption  of 
adiabatic  behavior  of  the  gas  in  the  bubble. 

2.  The  heating  due  to  viscous  dissipation  in  a  rapidly  compressed  drop  accounting  for  the  temperature 
dependence  of  the  liquid  viscosity. 

According  to  Field  (1992;  Field  et  al.  1992),  the  “hot  spots”  responsible  for  ignition  need  to  have 
dimensions  of  typically  0.1  -  10  /im,  duration  of  0.01  to  1  msec,  and  temperatures  greater  than  ~  700  K.  It 
is  useful  to  keep  in  mind  these  orders  of  magnitude  in  the  light  of  the  results  described  below. 


2  Bubble  Collapse 

One  of  the  scenarios  envisaged  for  the  observed  accidental  ignition  is  as  follows:  A  gas  bubble  can  remain 
trapped  by  several  mechanisms  in  the  course  of  propellant  transport  and  filling  operations.  When  this  bubble 
is  exposed  to  high  pressures  (e.g.  due  to  shock  waves  generated  by  water  hammer  effects  or  other  causes),  it 
collapses  and  the  nearly-adiabatic  heating  of  the  gas  can  cause  a  temperature  rise  of  the  surrounding  liquid 
sufficient  to  cause  ignition.  The  concrete  possibility  of -explosive  sensitization  due  to  gas  cavities  has  been 
demonstrated  for  solid  propellants  (see  e.g.  Heavens  and  Field  1974;  Chaudhri  and  Field  1974;  Field  et  al. 
1982;  Krishna  Moan  and  Field  1984).  For  liquid  propellants  the  experimental  evidence  is  less  clear  (Coley 
and  Field  1973;  Bourne  and  Field  1989,  1991,  1992).  While  sometimes  initiation  is  observed,  many  times  it 
is  not  and,  furthermore,  the  role  played  by  the  impacting  of  the  microjet  on  the  opposite  side  of  the  bubble 
is  unclear.  Theoretically,  the  situation  is  not  much  clearer.  While  some  authors  (e.g.  Andersen  and  Gillespie 
1980)  discount  the  role  of  gas  pocket  compression,  others  (Johansson  1958;  Randolph  and  Simpson  1976; 
Morrison  et  al.  1982;  Frey  1985)  reach  the  opposite  conclusion. 

In  all  of  the  previous  theoretical  analyses,  the  gas  behavior  in  the  collapsing  bubble  was  modelled  very 
crudely,  usually  by  assuming  adiabatic  behavior.  The  purpose  of  this  part  of  the  work  is  to  relax  this 
assumption  and  to  account  properly  for  the  thermal  behavior  of  the  gas.  We  shall  consider  two  cases, 
spherical  collapse  and  asymmetric  collapse  in  the  vicinity  of  a  plane  rigid  wall. 


2.1  Mathematical  Formulation 


We  assume  that  the  liquid  is  incompressible  and  inviscid  so  that  its  velocity  field  u  can  be  described  in  terms 
of  a  harmonic  potential  <f> 

u  =  V^,  V20  =  0.  (2.1.1) 


With  the  neglect  of  body 
Bernoulli  integral  is 


force  that  do  not  influence  the  process  of  present  concern  due  to  its  rapidity,  the 


9<j>  1  2  .  Pl_  _  P°°  (0 

dt  2 U  pl  Pl 


(2.1.2) 


where  pL  is  the  liquid  density,  Pl  the  pressure,  and  px  the  time-dependent  “ambient”  pressure  (the  rising  of 
which  will  cause  the  bubble  to  collapse).  The  kinematic  boundary  condition  on  the  bubble  surface  requires 
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that 


n  •  =  n  •  ,  (2.1.3) 

where  n  is  the  unit  normal  directed  into  the  liquid  and  x,  the  generic  surface  point.  With  the  neglect  of 
viscous  stresses  in  the  gas,  the  dynamic  boundary  condition  stipulates  that,  on  the  bubble  surface  S, 

p  —  Pl  —  &C  (2.1.4) 

where  p  is  the  gas  pressure,  a  the  surface  tension  coefficient,  and  C  the  local  surface  curvature  reckoned  as 
positive  when  the  centers  of  curvature  lie  in  the  gas.  In  particular,  this  relation  establishes  the  value  of  the 
internal  gas  pressure  pgo  necessary  for  a  spherical  bubble  with  radius  Ro  to  be  in  equilibrium  in  a  constant 
pressure  field  po: 

Pgo-Po  =  ^-  (2.1.5) 

In  our  treatment  of  the  gas  we  assume  that  the  Mach  number  of  the  gas  flow  is  small  so  that  the  gas 
pressure  p  is  approximately  uniform  in  space, 


p(x,t)  ~ p(t ). 


(2.1.6) 


This  assumptions  is  certainly  justified  during  the  greatest  part  of  the  collapse.  In  the  very  latest  stages  the 
bubble  wall  can  reach  velocities  of  several  thousand  m/s,  but  on  the  other  hand  the  gas  temperature  is  also 
quite  high  so  that  the  speed  of  sound  can  easily  be  greater  than  1,000  m/s.  Thus,  while  one  expects  (2.1.6) 
to  be  in  error  in  the  final  stages  of  the  collapse,  this  error  is  in  fact  smaller  than  might  appear  at  first  and 
is  also  significant  only  for  a  minute  fraction  of  the  total  process.  The  approximation  (2.1.6)  avoids  the  need 
to  consider  the  momentum  equation  in  the  gas. 

The  gas  continuity  equation  is 

^+pGV-v  =  0,  (2.1.7) 

where  pa  and  v  are  the  density  and  velocity  fields  in  the  gas  and  d/dt  =  d/dt  +  v  •  V  denotes  the  convective 
derivative,  and  the  gas  energy  equation  is 


(2.1.8) 


where  Cp  and  K  are  the  constant-pressure  specific  heat  and  thermal  conductivity  and  T  is  the  temperature 
field.  For  a  perfect  gas 

p  =  RpgT ,  (2.1.9) 

where  1Z  is  the  universal  gas  constant  divided  by  the  gas  molecular  weight,  so  that 


(2.1.10) 


Furthermore,  upon  combining  (2.1.7)  and  (2.1.8)  and  using  the  fact  that 


CppgT  =  - -p, 


where  7  is  the  ratio  of  specific  heats,  one  finds 

dp 
dt 


1  r  > 

7-1 

+  7pV  •  v  =  (7  -  1)V  •  (KVT) . 


Upon  using  the  assumed  spatial  uniformity  of  p,  this  equation  can  be  cast  in  the  form 

V  •  ^7pv  +  ^px  -  (7  -  1  )KVTj  =  0 , 


(2.1.11) 


(2.1.12) 


(2.1.13) 
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where  x  is  the  position  vector  and  the  dot  denotes  the  time  derivative.  This  equation  can  be  integrated  to 
find 

7 pv  +  -px  -  (7  -  \)I<VT  =  A ,  (2.1.14) 

where 

VA  =  0-  (2.1.15) 

Upon  solving  for  v  and  calculating  the  vorticity  field  we  find 

7pV  x  v  =  V  x  A +  (7  -  1)V  x  (ATVT).  (2.1.16) 

In  gases  the  thermal  conductivity  is  very  nearly  dependent  only  on  the  temperature,  so  that  the  last  term 
is  negligible.  This  relation  then  shows  that  VxvocVxA.  On  the  other  hand,  it  is  well  known  from 
Fluid  Mechanics  that,  in  a  motion  started  from  rest,  vorticity  can  only  arise  due  to  the  barotropic  term 
proportional  to  VpG  x  Vp  or  to  diffusion  from  the  boundaries  and  convection.  The  barotropic  term  vanishes 
here  due  to  (2.1.6)  and,  due  to  the  small  value  of  the  dynamic  viscosity  for  gases,  we  shall  neglect  boundary 
vorticity  diffusion.  As  a  consequence,  we  must  require  that  V  x  A  ~  0.  We  thus  may  write  A  =  with 
from  (2.1.15), 

V2$  =  0.  (2.1.17) 

Hence,  from  (2.1.14)  the  gas  velocity  field  can  be  expressed  as 

(7-l)A-Vr+V*-ipx]  •  (2.1.18) 

Another  consequence  of  (2.1.13)  is  an  equation  for  the  gas  pressure.  To  find  this  relation  we  integrate 
(2.1.13)  over  the  bubble  volume  V  and  use  the  divergence  theorem  to  find 

js  (iP*  +  ipx  -  (7  -  l)AVr)  •  n dS  =  0 ,  (2.1.19) 

where  S  is  the  bubble  surface.  From  this  relation  we  readily  obtain 

P  =  [(7  ~  1)Q  +  7 pV]  (2.1.20) 

where 

Q  =  -J  Kn  VTdS 

is  the  total  heat  conducted  into  the  liquid  at  the  bubble  wall. 

From  Eq.  (2.1.18)  and  the  kinematic  condition  we  find  a  boundary  condition  for  $  at 

n  •  •  n  +  (7  -  1)?„  +  |px,  •  n ,  (2.1.22) 

where 

?„  =  -ATn.Vr  (2.1.23) 

is  the  heat  flux  into  the  liquid  at  the  bubble  wall. 

One  last  equation  is  needed  to  complete  the  formulation,  and  this  is  taken  to  be  the  energy  equation 
(2.1.8)  specialized  to  the  case  of  a  perfect  gas: 

f  (f  +  Y  ' V!r)  -  *  =  V '  ■  (2-1-24) 


(2.1.21) 

the  bubble  wall: 
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it  al.  1988;  Prosperetti  1991).  (2.1.3)  we  have 

For  the  problem  in  the  liquid,  from  (2.1.1J  a  l  J  ^ 


<j>  =  -  —  R, 

v  r 


SKSSlS'--— : 

he^well-known  Rayleigh-Plesset  equation  f  ^  (,2,2) 


„  1  2cr"| 

RR  +  \R?  =  —  [P  -  p~®  -  "rJ 


«  J  ffp.ts  in  the  liquid  compressibility  are  readily  incorporated  (see 

F  v  r  a  -n  a 


(  r\  ..  3  ( .  R 

'-rjRR+i>\  ^ 


II  R  L 

R2  =  "7“  I  *  cl  cl  dt 


PL 


{,->  -W-^)-  (2M) 


l  Cr  I  l  J  u  -  .  .  , 

V  1  ,  „  fnr  S,  bounded  at  the  origin  is  $  - 

“d,«  asss  -  ~ lt“  “nsl“l  to  be  0  H““’ 

,e  gas  velocity  field  becomes,  from  (2.1.  ).  ^  ^  ,  (2,2.4) 


dT  1.1 

^TP[(^i)Av-3pi’ 

.  r  v  The  pressure  equation  (2.1.20)  gives  in  this  case 
iere  v  is  the  radial  component  of  v.  The  pressure  1 

3  r . dT\  _pl 


R 


b-^Tr 


I  r=H 


pon  substitution  into  (2.2.4)  we  have 
he  energy  equation  (2.1.24)  is 


7  p  fdT  dT' 
_.7—  —  I  — — (■  «-r 
T  -  1  T  \dt  dr 


p  “  r2  dr  V  5r  / 


(2.2.5) 


(2.2.6) 


(2.2.7) 


1  rp  \  fit  or  J  1  ' 

^  o  oiinw  the  thermal  conductivity 

rce  the  gas  temperature 

depend  on  temperature  and  rewrite  (^)  ^  } 


=  f  K{8)  d6  , 

J  Too 


,  rS  re  linon  T  has  been  assumed 
ere  T.  is  the  initial  temperature.  A  linear  relation  for  the  epen  ence  ^  ^  ^  ^  that 

hThXtraU, “symmetric  case  it  is  also  relatively  errsy  to  solve  the  energy  eq 

K  -  -  o  /  at,  \  (2.2.9) 


m.  r1  ksn  ^ 

-Qf+7?R  dr  r 2  dr  \  dr  J 
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where  Dl  is  the  liquid  thermal  diffusivity. 

An  efficient  numerical  model  for  the  solution  of  these  equations  has  been  described  in  Kamath  and 
Prosperetti  (1989)  and  Kamath  et  al.  (1993).  The  methd  is  based  on  a  coordinate  transformation  in  the  gas 


V  = 


m 


(2.2.10) 


and  in  the  liquid 


(2.2.11) 


and  a  spectral  representation  of  the  gas  and  liquid  temperature  fields.  In  (2.2.11)  Ro  is  the  initial  bubble 
radius,  u  a  characteristic  inverse  time,  and  £  a  numerical  constant  taken  to  be  about  0.05. 


2.3  Results  for  Spherically  Collapsing  Bubbles 

We  now  present  some  numerical  results  for  spherical  bubbles.  In  the  examples  that  follow,  unless  otherwise 
indicated,  the  physical  properties  of  air  are  used.  For  the  liquid  we  take  pL  =  1,452  kg/m3,  K  =  0.15  W/mK, 
Dl  =  4.49  xlO-8  m2/s-  These  values  are  representative  of  LGP  1845. 

The  equilibrium  pressure  is  po  =  1  bar  and  the  undisturbed  temperature  =  293  K. 

We  first  consider  Gaussian  pressure  pulses  with  the  form 


Poo  =  Po  +  A p  exp  , 


(2.3.1) 


where  A p  is  the  maximum  amplitude  and  tp  is  a  measure  of  the  width  of  the  pressure  pulse. 

Our  first  case  is  a  bubble  with  Ro  =  1  mm  for  A p  =  10  bars.  The  corresponding  value  of  the  Rayleigh 
collapse  time  tR  defined  by  _ 

tR  =  0.915Ro./^F  (2.3.2) 

V  AP 

is  29  ps.  Figure  1  shows  the  normalized  radius  R/Ro  versus  the  normalized  time  t/tR  for  tR/tp  =  0.1 
(dashed),  l(solid),  and  10  (dotted).  In  the  first  case  the  pressure  rise  is  slow  and  the  bubble  is  compressed 
quasi-statically  except  for  a  series  of  weak  nearly  free  oscillations  executed  about  the  instantaneous  pressure. 
When  the  pressure  pulse  is  much  shorter  than  the  collapse  time  (dotted  line),  on  the  other  hand,  the  bubble 
responds  with  essentially  free  oscillations  about  the  undisturbed  constant  pressure  level.  When  tR/tp  —  1, 
the  collapse  is  very  violent  with  a  volume  reduction  of  nearly  two  orders  of  magnitude. 

The  temperature  attained  by  the  gas  at  the  bubble  center  (squares)  and  by  the  liquid  at  the  bubble 
surface  (circles)  at  the  end  of  the  first  collapse  are  shown  in  Fig.  2  for  tR/tp  =  0.1,  0.2,  0.5,  1,  2,  5,  and  10. 
The  circles  in  this  figure  are  the  adiabatic  prediction  for  the  gas  temperature 


(2.3.3) 


It  can  be  seen  that  the  heating  is  maximum  when  the  pressure  time  scale  tp  is  close  to  the  bubble  collapse 
time  tR.  In  none  of  these  cases,  however,  the  liquid  heating  appears  to  be  sufficient  to  cause  ignition. 

The  center  gas  temperature  and  the  liquid  surface  temperature  for  the  three  cases  of  Fig.  1  are  shown 
as  a  function  of  t/tR  in  Figs.  3  and  4  respectively  (the  figures  are  drawn  to  the  same  scale).  As  before  the 
dashed  curves  are  for  tR/tp  =  0.1,  the  solid  ones  for  tR/tp  =  1,  and  the  dotted  ones  for  tR/tp  =  10.  The 
marked  differences  caused  by  the  “detuning”  of  the  rate  of  compression  with  respect  to  the  collapse  time 
are  very  clear  from  these  figures  also.  Figure  5  shows  the  corresponding  wall  heat  flux.  The  duration  of  the 
temperature  spike  for  tR/tp  =  1  is  less  than  2  ps  and  therefore  shorter  than  the  lower  bound  on  the  duration 
required  for  ignition  quoted  in  section  1.  The  maximum  temperature  reached  by  the  liquid,  on  the  other 
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hand  (Fig.  4)  seems  to  be  too  small  to  cause  ignition.  Figure  6  is  the  amount  of  heat  Q  delivered  to  the 
liquid  as  a  function  of  time: 

dt .  (2.3.4) 

r=R(t) 

Since  these  results  show  that  the  largest  temperature  rises  occur  for  tR~tp,  we  now  consider  the  effect 
of  the  bubble  radius  when  tR  =  tp  for  the  same  A p  =  10  bars  in  Figs.  7  and  8.  These  figures  show  the  core 
gas  (dashed)  and  the  bubble  surface  (solid)  temperature  for  Ro  =  100  pm  and  10  mm  respectively.  The 
corresponding  values  of  tR  are  2.9  ps  and  0.29  ms.  The  temperature  responses  are  very  similar  when  plotted 
as  functions  of  the  normalized  time  t/tR  as  here.  The  fifst  compression  gives  a  high  temperature  in  the  gas, 
but  the  rise  in  the  liquid  is  relatively  modest.  The  subsequent  gas  temperature  oscillations  due  to  the  free 
oscillations  executed  after  the  passage  of  the  pulse  (see  Fig.  1)  hardly  afFect  the  liquid  temperature. 

When  the  magnitude  of  the  pressure  pulse  is  larger,  A p  =  100  bars,  as  in  Figs.  9  to  11,  the  general 
behavior  is  also  little  dependent  on  the  value  of  the  radius  on  a  t/tR  scale.  These  figures  refer  to  Ro  =  0.1,  1, 
and  10  mm.  The  corresponding  values  of  tR  are  0.92,  9.2,  and  92  ps.  The  liquid  temperature  accompanying 
the  large  gas  temperature  spike  reaches  about  1,500  K,  but  the  duration  of  this  extreme  temperature  seems 
to  be  too  short  to  cause  ignition  even  for  the  10  mm  bubbje  for  which  it  is  of  the  order  of  5  ps.  Furthermore,  it 
appears  likely  that  the  liquid  temperature  would  be  somewhat  decreased  by  vaporization  that  is  not  included 
in  the  present  model. 

Another  situation  of  interest  is  that  in  which  an  overpressure  of  long  duration  is  applied  to  the  bubble: 

PooOO  =  po  +  ApH(t)  (2.3.5) 

where  H(t)  is  Heaviside’s  step  function.  This  simulates  the  effect  of  a  compression  wave  with  a  rise  time 
faster  than  the  Rayleigh  time  and  a  duration  much  longer  than  tR.  For  a  slower  rise  time  the  bubble  would 
compress  quasi-statically  as  shown  before  in  Figs.  1,  3,  and  4  and  the  temperature  increases  would  be  small. 

Figures  12  to  15  are  for  A p  =  10  bars.  Figure  12  shows  the  bubble  radius  versus  time  for  Ro  =  1  mm.  The 
pressure  step  excites  oscillations  that  are  slowly  damped  by  the  dissipative  processes  affecting  the  bubble 
motion.  Asymptotically,  the  bubble  radius  tends  to  a  new  equilibrium  value  consistent  with  the  pressure 
po  +  Ap,  see  (2.1.5).  This  behavior  is  typical  for  all  the  cases  studied.  Figures  13  to  15  show  the  temperature 
response  of  the  gas  and  of  the  liquid  for  Ro  =  0.1,  1  and  10  mm  respectively.  The  liquid  temperature  rises 
after  each  collapse  due  to  the  gas  temperature  spike,  but  the  integrated  effect  remains  small. 

The  radius-versus-time  curve  for  Ro  =  1  mm  and  a  larger  pressure  step,  Ap  =  100  bars,  is  shown  in  Fig. 
16.  The  gas  and  liquid  temperature  behavior  under  this  stronger  excitation  are  shown  in  Figs.  17  to  19 
again  for  Ro  =  0.1,  1  and  10  mm.  The  liquid  temperature  is  somewhat  greater  than  before,  but  still  appears 
to  remain  below  the  level  required  for  ignition  except  for  an  initial  spike  with  a  duration  of  less  than  10  ps. 
It  is  doubtful  that  such  marginal  conditions  would  be  sufficient  to  trigger  ignition. 

We  have  also  conducted  a  few  simulations  with  a  gas  adiabatic  index  y  =  5/3  (appropriate  for  monatomic 
gases)  and  different  values  of  the  liquid  thermal  conductivity.  We  give  an  example  in  Fig.  20  with  K  = 
0.075  W/mK  (i.e.,  half  the  value  used  in  the  previous  simulations)  y  =  5/3,  Ro  =  10  mm,  for  a  step  pressure 
increase  Ap  =  100.  This  figure  should  therefore  be  compared  with  Fig.  19.  Due  to  the  greater  value  of  y 
the  gas  heating  is  much  larger  than  before  and  is  not  shown  in  order  to  maintain  readability  of  the  liquid 
temperature  information.  The  liquid  gets  much  hotter  than  before  during  the  initial  spike  which  in  this  case 
is  somewhat  longer  than  before  (about  20  instead  of  10  pm).  Although,  on  the  basis  of  the  estimate  given 
in  the  Introduction,  ignition  would  appear  to  be  possible,  conditions  seem  to  be  too  close  to  marginal  to  be 
able  to  make  a  definite  prediction. 

2.4  Non-Spherical  Collapse  -  Full  Model 

In  the  case  of  a  non-spherical  collapse,  no  analytical  progress  is  possible  and  the  mathematical  model  outlined 
before  must  be  solved  by  fully  numerical  means. 
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surface  and  using  (2.1.4)  to  eliminate  the  liquid  pressure  we  find 


=  d±  +  un^  =  ul-W- 

nr  -  dt  +Un8n  "  2 


p  —  cC  —  Poo  (0 


(2.4.1) 


adapted  to  the  present  axisymmetnc  situation: 

4>x,t)  =  Js(Gun-H<t>)ds'.  (2-4-2) 


Here  the  integration  is  along  the  trace  £  of  the  bubbte  surface  on  * 

H  are  related  to  elliptic  mtegrals  and  depend  on  the  field  point  *  andm,  ^  ^  ^  fte  ^ 

fo-  l  ®  ven  ^  inown  ire.  integrating  (2.4.!). 

The  numerical  procedure  used  to  execute  a  time  step  can  be  summarized  as  follow  . 

1.  The  potential  on  the  bubble  surface  is  updated  by  the  explicit  one-step  Euler  method  using  (2.4.1) 

2.  The  velocity  normal  to  the  surface  is  then  evaluated  from  (2.4.2) 

3.  With  this  velocity  field  the  bubble  surface  is  displaced  using  (2.1.3),  again  with  the  explicit  one-step 
Euler  method 

4.  An  orthogonal  grid  is  generated  inside  the  bubble;  an  example  is  shown  in  Fig.  21. 

5.  The  Laplace  (2.1.17)  and  energy  (2.1.24)  equations  are  discrete 

finite  differences  and  the  resulting  system,  including  the  pressure  equal, on  (2.1.20),  ,s  solved  ,terat,ve,y. 

Unfortunately,  this  numerical  scheme  suffers  from  a  weakness  in  that,  as  is *  C'f °  the^aspect  ratio 

retrtht^ 

•SSHsSSSis: 

fllustration  of  this  feature  of  the  temperature  d.tnbutmn  ,s  proved  “Jf  ”  temperature  is 

near  the  boundary  to  formulate  an  approximate  calculat.on 

scheme  that  we  now  describe. 
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2.5  Non-Spherical  Collapse  -  Boundary  Layer  Approximation 

The  very  abrupt  change  of  the  temperature  field  near  the  bubble  boundary  shown  in  Fig.  23,  coupled  with 
the  smooth  variation  of  the  field  $,  suggests  to  approximate  the  energy  equation  (2.1.24)  in  a  boundary 

lay  To This^end ,  we  introduce  a  new  coordinate  system  in  which  the  bubble  surface  corresponds  to  the  line 
x  _  o  while  the  y  =  constant  lines  are  normal  to  the  surface.  In  the  immediate  vicinity  of  the  surface 
these  coordinate  lines  are  very  nearly  orthogonal.  In  this  system,  the  energy  equation  (2.1.24)  is,  omitting 
conduction  in  the  tangential  direction, 


T7  +  (®*  ~  -  vy^%  =  ~Di>  +  D W ' 


dt 


dy 


(2.5.1) 


Here  r  is  the  auxiliary  temperature  variable  defined  in  (2.2.8)  and  D  -  K/(Cppc)  is  the  gas  thermal 
diffusivity.  Furthermore,  the  term  Vx  is  the  grid  velocity  in  the  tangential  direction  and  i>y0  is  the  bubble 

surface  normal  velocity.  From  (2.1.18)  we  have 


P  ,  7~ 
Vy°  ~  3j PV  7 P 


y=0j 


<9$  5$ 

+  dy  dy 


y=0 


(2.5.2) 


In  view  of  the  smooth  variation  of  $  in  the  vicinity  of  the  bubble  wall  (see  Fig.  22)  the  last  two  terms  can 
be  simplified  by  means  of  a  Taylor  series  expansion: 


cM> 

dy 


d2$ 

W  y=0V 

~\7h*+Xdy  +  Ry  dy)  V' 


(2.5.3) 


where  Rx  and  Ry  are  the  radii  of  curvature  of  the  surfaces  x  =  const,  and  y  -const  calculated  from  the 
computed  surface  shape.  The  step  from  the  first  to  the  second  line  is  due  to  the  fact  that  V  9  -  0. 

The  numerical  solution  procedure  is  the  same  as  in  the  previous  section,  except  for  the  fact  that  no  gri 
is  generated  numerically  and  only  an  algebraic  stretching  is  introduce  to  increase  the  node  density  near  the 

This  boundary  layer  approximation  has  been  compared  with  the  full  model  results  for  the  spherical  case 
and  for  the  asymmetric  case  up  to  the  point  of  jet  formation  with  excellent  results.  A  typical  comparison 
for  the  spherical  case  is  shown  in  Fig.  24  where  the  dashed  line  is  the  gas  center  temperature  according  o 
the  full  model  and  the  solid  line  the  boundary  layer  approximation.  A  comparison  for  the  asymmetric  case 
is  shown  in  Fig.  25  that  shows  the  heat  fluxes  at  the  jet  tip  predicted  by  the  full  model  (dashed  line)  and 
by  the  boundary  layer  approximation  (solid  line)  for  the  case  of  Figs.  21  to  23. 


2.6  Non-Spherical  Collapse  —  Some  Results 


Figure  26  shows  successive  bubble  configurations  for  collapse  induced  by  a  sudden  overpressure  Eq.  (2.3.5), 
with  A p  =  10  bars,  Po  =  1  bar.  Here  the  initial  bubble  radius  is  Ro  =  1  mm  and  the  rigid  wall  is  initially 
1  2  mm  below  the  bubble  center.  Successive  bubble  shapes  for  the  excitation  by  the  pressure  pulse  (2.3.1) 
with  tR/tv  =  4  are  shown  in  Fig.  27.  The  peculiar  shjfre  of  the  free  surface  in  this  case  is  due  to  the  fact 
that,  once  formed,  the  jet  continues  its  motion  toward  the  opposite  wall  even  after  the  rest  of  the  cavity  as 

started  to  rebound.  _  .  .  in 

Figure  28  shows  the  local  heat  flux  along  the  bubble  surface  for  the  step  pressure  increase  A p  =  10  bars 

as  computed  from  the  boundary-layer  model  as  a  function  of  time.  The  abscissa  in  this  figure  denotes  the 
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node  number.  Sixty  equispaced  nodes  were  used.  The  heat  flux  is  fairly  uniform  over  the  surface  until  the  jet 
forms.  At  this  point  the  region  near  the  jet  tip  starts  absorbing  heat  considerably  faster  than  other  regions. 
Just  before  the  jet  strikes  the  other  side  of  the  bubble,  the  heat  flux  at  the  tip  is  calculated  to  reach  more 
than  300  MW/m2! 

In  this  case  we  were  unable  to  calculate  directly  the  temperature  in  the  liquid  to  find  the  temperature 
at  the  bubble  surface.  Hence,  we  use  the  following  approximation.  At  the  jet  tip,  the  flow  is  only  in  the 
direction  normal  to  the  axis,  and  in  the  immediate  vicinity  of  the  tip  (which  is  the  only  region  where  there 
are  appreciable  temperature  gradients)  is  very  nearly  solid  body.  In  a  frame  of  reference  moving  with  the 
tip,  therefore,  the  liquid  energy  equation  is  therefore  essentially  the  pure  conduction  equation.  Furthermore, 
conduction  in  the  direction  tangent  to  the  surface  is  very  likely  much  slower  than  in  the  thin  boundary  layer 
normal  to  the  surface.  These  considerations  satisfy  solving,  on  the  tip  axis,  the  equation 

=  (2.6.1) 
m  L  dyi  '  '  ’ 

subject  to  continuity  of  the  normal  heat  fluxes.  Even  aside  from  the  previously  described  approximations  this 
procedure  is  not  strictly  correct  as  the  gas  energy  equation  is  solved  assuming  that  the  liquid  temperature 
is  undisturbed.  However  the  gas  heating  is  so  large  that  it  is  unlikely  that  accounting  for  the  liquid  surface 
heating  would  affect  the  gas  temperature  field  very  much. 

In  this  way  we  can  estimate  the  liquid  temperature  at  the  jet  tip  as  a  function  of  time.  Results  for  the 
case  of  Fig.  26  are  shown  in  Fig.  29  (solid  line).  In  this  figure  the  dashed  line  is  the  liquid  temperature  fr 
the  corresponding  spherical  case. 

2.7  Conclusions  on  Bubble  Collapse 

We  have  presented  results  for  the  gas  and  liquid  temperatures  attained  during  the  collapse  of  gas  bubbles 
exposed  to  pressure  pulses  and  steps.  The  main  conclusions  are  as  follows: 

1.  The  ratio  between  the  characteristic  time  for  the  collapse  (the  Rayleigh  time)  Ir  and  the  duration  of 
the  pressure  pulse  tp  has  emerged  as  a  very  significant  parameter.  When  this  ratio  is  close  to  1  the 
gas  temperature  reaches  very  high  values.  When  the  pressure  pulse  amplitude  is  100  bars  and  tp/tR 
=  1  our  computations  give  peak  liquid  temperatures  of  the  order  of  1500  K  for  bubble  radii  between 
0.1  and  10  mm.  The  duration  of  these  large  temperatures  is  however  very  short,  most  likely  too  short 
to  result  in  ignition.  Furthermore,  it  may  be  assumed  that,  under  these  conditions,  phase  change  will 
occur  at  the  bubble  wall  with  a  corresponding  decrease  in  temperature. 

2.  We  have  also  studied  the  response  of  a  spherical  bubble  to  a  pressure  step.  In  this  case  the  bubble  ex¬ 
ecutes  a  series  of  free  oscillations  of  decreasing  amplitude.  As  before,  we  find  large  liquid  temperatures 
at  the  bubble  surface,  especially  during  the  first  collapse  that  is  the  most  violent  one.  These  tem¬ 
peratures  are  the  larger  the  smaller  the  liquid  thermal  conductivity.  Also,  other  things  being  equal,  a 
monatomic  gas  gets  hotter  than  a  diatomic  one  due  to  the  higher  adiabatic  index,  with  a  corresponding 
greater  liquid  temperature.  On  the  basis  of  these  previous  considerations,  one  may  perhaps  postulate 
a  new  mechanism  for  ignition:  upon  the  first  collapse,  the  bubble  surface  heats  up  considerably 
and  part  of  the  liquid  vaporizes  into  the  bubble  core.  Upon  subsequent  collapses  this  vapor  is  exposed 
to  very  large  temperatures,  much  larger  than  the  liquid  at  the  bubble  surface,  and  ignites.  Whether 
this  mechanism  is  possible  cannot  be  decided  on  the  basis  of  the  current  model  that  does  not  include 
phase  change  processes.  An  extension  in  this  direction  should  however  be  feasible. 

3.  All  of  the  previous  conclusions  have  been  reached  assuming  that  the  bubble  maintains  a  spherical  sym¬ 
metry.  This  is  actually  rather  unlikely  in  view  of  the  well-known  instability  of  the  spherical  shape.  For 
this  reason,  and  also  to  simulate  the  effect  of  nearby  boundaries,  we  have  formulated  a  model  capable 
of  handling  axisymmetric  collapses.  This  model  is  considerably  more  complex  than  the  symmetric  one 
and  only  a  limited  exploration  of  its  predictions  has  been  possible  whithin  the  time  span  of  this  grant. 
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Nevertheless,  as  shown  in  Fig.  29,  we  find  that  the  liquid  microjet  that  penetrates  the  bubble  as  a 
result  of  the  instability  is  exposed  to  much  higher  temperatures  than  in  the  spherical  case,  although 
again  the  duration  of  these  high  temperatures  is  very  short.  Experimental  evidence  does  indicate  that 
ignition,  when  it  occurs,  seems  to  be  related  (at  least  in  time)  to  the  impact  of  the  microjet  on  the 
other  side  of  the  bubble  wall.  It  is  not  clear  however  whether  this  is  the  result  of  the  microjet  heating 
up,  or  of  the  fact  that  the  collapse  continues  even  after  the  impact  of  the  microjet,  and  that  the  maxi¬ 
mum  compression  (and  therefore  the  strongest  heating)  immediately  follows  the  microjet  impact.  This 
circumstance  is  really  accidental,  in  the  sense  that  it  is  only  due  to  the  fact  that  the  microjet  takes  a 
very  long  time  to  develop. 

4.  All  of  the  results  that  we  have  found  indicate  that  ignition  is  in  principle  possible,  although  marginally 
so.  This  finding  may  explain  why  it  is  difficult  to  deduce  definite  conclusions  from  the  available 
experimental  evidence  with  its  apparently  inconsistent  results. 
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21.  Example  of  a  numerically  generated  orthogonal  grid  for  the  solution  of  the  Laplace  and  energy 
iion  in  the  gas.  Here  the  initial  bubble  has  a  radius  R o  =  1  mm,  with  a  rigid  wall  1.2  mm  below  the 
I  bubble  center.  The  collapse  is  triggered  by  a  pressure  step  A p  =  10  bars.  The  grid  is  shown  at  the 
t  =  26.64  (is.  The  Rayleigh  time  is  in  —  29  /js. 


le  left 


Temperature  Profile 


Fig.  23.  The  gas  temperature  along  the  axis  of  symmetry  of  the  bubble  (from  the  closest  to  the  farthest 
point  from  the  wall)  for  the  case  of  the  previous  figure. 
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Estimated  jet  tip  liquid  temperature  as  a  function  of  time  for  the  case  of  Figs.  21  to  23  (solid 
line).  The  dashed  line  is  the  liquid  temperature  for  a  spherical  collapse. 


3  Viscous  heating  of  liquid  layers  under  intense  shear 

The  second  mechanism  that  might  lead  to  the  appearance  of  a  “hot  spot”  in  the  liquid  monopropellant  that 
we  consider  is  viscous  heating,  with  particular  regard  to  the  flow  occurring  during  the  squeezing  of  a  drop 
of  propellant  between  two  colliding  parallel  plane  surfaces.  This  particular  configuration  is  suggested  by  the 
extensive  experimental  work  carried  out  by  Field  and  co-workers  in  their  falling  weight  apparatus  (Heavens 
and  Field  1974;  Swallowe  and  Field  1982;  Krishna  Moan  and  Field  1984;  Krishna  Moan  et  al.  1984;  Field  et 
al.  1982;  Field  et  al  1992;  Bourne  and  Walley  1993).  Although  this  situation  has  been  studied  before  (Eirich 
and  Tabor  1947;  Field  et  al.  1982;  Bourne  and  Walley  1993),  the  effects  of  the  temperature  dependence 
of  viscosity  have  been  either  disregarded  or  treated  only  approximately.  Yet,  this  feature  of  the  problem 
is  exceedingly  important  as  shown  by  the  estimates  of  maximum  temperature  rise  given  by  Heavens  for 
nitroglycerine  (Field  et  al.  1982)  which  drop  from  about  650  K  to  100  K  when  the  temperature  dependence 
of  viscosity  is  approximately  allowed  for. 

In  the  flow  considered,  the  maximum  temperature  rise  occurs  in  thin  regions  parallel  to  the  colliding 
surfaces.  In  these  regions  the  viscosity  is  minimum  and  the  shear  is  large.  There  is  therefore  a  resemblance 
with  the  occurrence  of  shear  bands  in  solid  materials,  which  are  known  to  be  one  of  the  processes  leading 
to  ignition  (Winter  and  Field  1975;  Field  et  al.  1982;  Field  1992).  This  feature  of  the  flow  emerges  very 
clearly  from  our  numerical  calculations,  but  is  shown  to  decrease,  rather  than  increase,  energy  dissipation. 
A  rather  non-conventional  mechanism  leading  to  another  analogue  of  a  shear  band  that  we  investigate  is  the 
allowance  of  a  small  amount  of  slip  at  the  solid-liquid  interface.  The  process  is  very  sensitive  to  this  effect, 
which  is  however  limited  by  the  expected  smallness  of  the  slip  coefficient.  This  conclusion  is  however  only 
tentative  as  little  quantitative  information  is  available  on  this  quantity. 


3.1  Mathematical  Model 


In  the  experimental  situation  to  which  reference  was  made  earlier,  the  portion  of  the  liquid  region  not 
bounded  by  solid  surfaces  is  in  contact  with  air.  During  the  flow,  therefore,  the  liquid-solid-air  contact  line 
moves  along  the  solid  surfaces.  The  problems  associated  with  a  realistic  modeling  of  this  motion  are  well 
known  (see  e.g.  Dussan  V.  1979).  Since,  in  this  study,  we  wish  to  focus  specifically  on  viscous  heating,  we 
have  decided  to  avoid  any  contact  line  problem  by  supposing  that  the  two  colliding  plates  are  immersed  in 
liquid  rather  than  air  as  shown  in  Fig.  1.  This  figure  also  shows  the  frame  of  reference  centered  on  the  axis 
of  symmetry  at  the  midpoint  of  the  gap  with  instantaneous  thickness  2 h(t).  The  radius  of  the  plates  is  R. 

The  flow  is  assumed  to  be  axi-symmetric  and,  since  the  distance  between  the  plates  is  much  smaller  than 
R,  we  make  a  boundary  layer  approximation  assuming  that  derivatives  in  the  radial  direction  r  are  small 
with  respect  to  those  in  the  axial  direction  normal  to  the  plates  z.  Neglecting  further  the  compressibility  of 
the  liquid,  the  equations  to  be  solved  are  the  equation  of  continuity 


1  dru  dw 
r  dr  dz 


(3.2.1) 


of  radial  momentum 


and  of  energy 


f  du  du  du\  dp  d  f  du\ 
'1ST  dT  dT\  ,  /SuN’ 

'’cHaT  +  “a7  +  u’e7j  =  a?  +  "  (s7j 


(3.2.2) 

(3.2.3) 


Here  the  symbols  have  their  customary  meaning:  u  and  w  are  the  velocity  components  in  the  radial  and 
axial  direction,  p  is  the  pressure,  T  is  the  temperature,  p  is  the  density,  Cp  is  the  specific  heat,  and  p  is  the 
viscosity  for  which  we  assume  a  temperature  dependence  of  the  form 


T 

p=  Poo  exp  T°Ti 


(3.2.4) 
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In  keeping  with  the  boundary  layer  approximation  the  pressure  is  taken  to  be  uniform  in  the  z  direction 
and  to  only  depend  on  the  distance  from  the  axis  of  symmetry,  p  =  p(r,t). 

The  flow  region  considered  is  a  cylinder  of  radius  R  bounded  by  the  plates.  The  lateral  surface  of  this 
cylinder  is  an  ideal  surface  immersed  in  the  fluid  on  which  we  prescribe  a  constant-pressure  condition, 

p  =  0  r  =  R.  (3.2.5) 


In  view  of  the  parabolic  nature  of  the  equations,  no  further  conditions  on  velocity  or  pressure  are  necessary 
at  r  =  R. 

The  normal  velocity  relative  to  the  plates  must  vanish, 


w(r,z  =  ±h,t)  =  TV(t)  on  z  —  ±A(t) , 


(3.2.6) 


where  V(t)  >0  is  the  upward  velocity  of  the  lower  plate  or,  equivalently,  the  downward  velocity  of  the  upper 
plate.  For  the  tangential  velocity,  we  write 

u  =  on  z  —  ±h(t) ,  (3.2.7) 

oz 

where  A  is  a  constant  with  dimensions  of  a  length.  Upon  taking  A  =  0,  we  recover  the  standard  no-slip 
boundary  condition. 

In  the  course  of  the  experiment  that  we  model,  the  plates  are  driven  together  by  the  impact  of  a  weight. 
When  the  mass  M  of  this  weight  is  very  large,  one  may  assume  that  the  velocity  V  is  a  constant.  However, 
for  smaller  masses,  V(t)  must  be  determined  accounting  for  the  normal  force  exerted  on  the  plates  by  the 
squeezing  of  the  liquid.  Thus,  in  this  case,  we  write 


(3.2.8) 


where  the  integration  is  over  the  surface  of  the  plate  from  0  to  R.  The  contribution  from  the  normal  viscous 
stress  is  identically  zero  in  the  no-slip  case  as  a  consequence  of  (3.2.1)  and  very  small  in  the  case  with  slip 
and  has  been  neglected.  The  factor  of  2  in  the  left-hand  side  accounts  for  the  fact  that  the  plate  velocity  in 
the  laboratory  frame  where  the  lower  plate  is  at  rest  is  twice  the  plate  velocity  in  the  present  frame. 

For  the  temperature  field  we  consider  two  conditions  on  the  solid  surfaces: 


i  Isothermal  wall 


T  =  T0. 


(3.2.9) 


ii  Adiabatic  wall 


We  now  introduce  dimensionless  variables,  denoted  by  an  asterisk,  as  follows: 


(3.2.10) 


r=Rr.  z  =  h(t)z .  h(t)  =  h0h* (t) ,  (3.2.11) 

Vo 

u  =  Vbu*  w  =  —. Vow. ,  (3.2.12) 

K 

p  =  pV02p . ,  T  =  ToT. ,  p  =  pop,  •  (3.2.13) 

Here  h0  =  /i(0)  is  one  half  of  the  initial  plate  separation  and  p0  is  the  liquid  viscosity  at  the  initial  ambient 
temperature.  The  definition  of  z.  in  (3.2.11)  has  the  advantage  of  mapping  the  time  dependent  domain 
—h(t)  <  z  <  h(t)  onto  the  fixed  domain  — 1  <  z*  <  1.  Vo  is  a  characteristic  value  of  the  velocity  V  with 
which  the  plates  approach  each  other  and  the  scaling  (3.2.12)  for  w  reflects  the  boundary-layer  like  nature 
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of  the  flow.  In  terms  of  these  new  variables,  and  dropping  the  asterisks  for  simplicity,  Eqs.  (3.2.1)  to  (3.2.3) 
become 

7^+^!r  =  °-  <«14> 
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Here  h!  —  dh/dt, 


p  =  exp 


7?(1-T) 


(3.2.15) 

(3.2.16) 

(3.2.17) 


(r — rt*)(i  —  Tj )  ’ 

and  the  aspect  ratio  A,  Reynolds  number  Re,  Prandtl  number  Pr,  and  Eckert  number  Ec  are  defined  by 


.  ho 

A  =  H' 


Re  = 


pVoho 
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Pr  = 


_  CpflQ 
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Ec  = 


Vo2 

CPT0 


(3.2.18) 


An  expression  for  the  pressure  gradient  can  be  found  by  integrating  Eq.  (3.2.15)  across  the  half-width  of 
the  gap  0  <  z  <  1: 


jdp  _  1  p  du 


dr  Re  h(t)2  dz 


d  f1  ,  A  fh'V  A  f1  dru2  ,  /ooin, 

,  -  m  l udz +  2  (tt)  r  -  7 1  ■  (3'219) 


Because  of  the  symmetry  of  the  flow  about  the  midplane  of  the  gap,  it  is  sufficient  to  solve  Eqs.  (3.2.14) 
to  (3.2.16)  in  the  range  0  <  z  <  1.  Appropriate  symmetry  conditions  are  imposed  on  z  —  0.  At  z  =  1  the 
velocity  condition  (3.2.7)  becomes 


\  du  tv,  \  X 

u  =  A,  — —  with  A,  =  — 

oz  ho 


and  the  dimensionless  thermal  boundary  conditions  are 
i  Isothermal  wall 


T  =  1. 


ii  Adiabatic  wall 


i  =  EcPr: 


A pfduy 
h  \dz ) 

Since  2V^  =  —dh/dt,  from  Eq.  (3.2.8)  we  have 

d2h  m  f1.  2\&Pj 


(3.2.20) 

(3.2.21) 

(3.2.22) 

(3.2.23) 


where  m  =  pirR2h0,  is  the  initial  mass  of  the  liquid  in  the  volume  of  interest.  In  deriving  this  equation  we 
have  expressed  p  in  (3.2.8)  as  the  integral  of  its  derivative  and  interchanged  the  order  of  integration. 
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3.2  Numerical  Method 

The  system  of  equations  posed  in  the  previous  section  is  solved  by  finite  differences  in  the  region  0  <  r,  z  <  1. 
In  order  to  cluster  the  grid  points  near  the  solid  surface  where  gradients  are  large,  we  use  a  uniform  grid  in 
the  auxiliary  variable  £  =  tan  (j ez)  /  tan  (^e),  that  also  ranges  between  0  and  1. 

We  discretize  the  momentum  and  energy  equations  (3.2.15)  and  (3.2.16)  fully  implicitly  using  upwinding 
for  the  convective  terms.  Upwinding  is  also  used  in  the  last  term  of  the  pressure  equation  (3.2.19).  The 
resulting  nonlinear  system  is  solved  iteratively  in  the  following  fashion.  First,  with  the  most  recently  updated 
value  of  the  temperature  the  viscosity  is  calculated  and  the  velocity  field  ti  updated  from  (3.2.15).  The 
pressure  term  in  the  equation  is  expressed  by  means  of  (3.2.19)  with  the  integrals  effected  by  the  trapezoidal 
rule.  Next  w  is  calculated  from  the  continuity  equation  (3.2.14)  and  the  temperature  field  updated.  The 
process  is  arrested  when  the  maximum  relative  correction  to  the  temperature  field  is  less  than  10-8. 

When  Eq.  (3.2.23)  for  V  is  used,  it  is  integrated  explicitly  by  the  Euler  one-step  method. 

To  simulate  the  impulsive  nature  of  the  flow,  initially  we  assume  that  the  radial  velocity  field  u  is  uniform 
and  the  vertical  velocity  w  is  linear  in  the  vertical  coordinate.  These  fields  are  such  that  the  continuity 
equation  and  the  normal  velocity  condition  (3.2.6)  are  satisfied  for  a  plate  velocity  V(0).  This  state  of 
motion  violates  of  course  the  no-slip  condition  but  is  an  accepted  procedure  for  impulsively  started  flows 
(see  e.g.  Gresho  1991). 

The  time  step  is  set  proportional  to  the  gap  width  to  resolve  the  extremely  large  velocity  and  temperature 
rates  of  change  in  the  final  stages  of  the  process. 

3.3  Results 

The  numerical  results  that  follow  have  been  obtained  for  the  case  of  LGP  1845  for  which  p  =  1.452  g/cm3, 
K  =  0.15  W/mK,  Cp  =  2,300  J/kgK,  Pr  =  175,  Ta  =  524  K,  Tb  =  167  K,  p^  =  0.169  cP.  We  also  take  the 
plate  radius  R  —  10  mm,  half  the  initial  plate  separation  ho  =  1  mm,  and  the  initial  temperature  T0  =  20 
°C,  with  a  corresponding  viscosity  po  =  10  cP.  With  an  initial  velocity  V(0)  =  1  m/s,  we  thus  have 

A  =  0.1,  Re  =  134.1,  Ec  =  1.4  x  10"6 .  (3.4.1) 

We  first  study  the  case  in  which  the  velocity  of  the  plates  is  prescribed  and  their  surface  is  adiabatic. 

In  view  of  the  cylindrical  geometry,  radial  velocities  get  larger  with  increasing  distance  from  the  axis  and 
we  therefore  focus  on  the  velocity  and  temperature  distributions  at  r  =  R.  In  Fig.  2  we  show  the  distribution 
of  u  across  the  gap  0  <  z/h(t)  <  1  at  different  times  for  V(0)  =  1  m/s.  The  calculation  is  stopped  when 
the  gap  thickness  reaches  10  pm,  i.e.  0.5%  of  the  initial  value.  Although  not  very  clear  from  the  figure, 
initially  the  velocity  distribution  is  very  nearly  uniform  except  for  a  thin  boundary  layer  near  the  plates. 
With  time,  viscous  diffusion  promotes  the  appearance  of  a  profile  with  a  vague  similarity  to  the  parabolic 
distribution  that  would  be  given  by  a  lubrication  approximation.  Since  here  the  Reynolds  number  for  the 
radial  flow  is  of  the  order  of  104  and  the  flow  highly  unsteady,  the  physics  underlying  this  distribution  is  of 
course  quite  different  from  that  applicable  in  the  lubrication  case.  Still  later,  due  to  the  increasing  energy 
dissipation  near  the  wall,  viscosity  decreases  and  the  velocity  becomes  nearly  uniform  in  the  central  region. 
As  a  consequence,  the  shear  becomes  localized  near  the  wall,  where  the  temperature  grows  at  a  faster  rate. 
Viscosity  decreases  further,  and  a  non-monotonic  velocity  profile  develops  exhibiting  features  analogous  to 
those  of  shear  band  in  a  solid.  In  the  example  studied  here,  this  feature  arises  only  very  late  in  the  process 
when  the  plate  spacing  is  of  the  order  of  0.5%  of  its  initial  value. 

A  corresponding  sequence  of  temperature  distributions  across  the  gap  at  r  =  R  is  shown  in  Fig.  3.  In 
view  of  the  large  value  of  the  Prandtl  number  the  temperature  rise  occurs  in  a  region  that  remains  relatively 
thin  compared  with  the  gap  width.  The  temperature  increase  becomes  substantial  only  in  the  very  last 
stages  of  the  process.  This  point  is  clearer  from  Fig.  4  which  shows  the  maximum  temperature  in  the  system 
(that  occurs  in  this  case  on  the  plates  at  r  =  R)  as  a  function  of  the  dimensionless  gap  width  h(t)/h( 0).  The 
corresponding  history  for  the  radial  velocity  at  r  =  R,  z  =  0  is  shown  in  Fig.  5.  In  these  two  figures  the  solid 
lines  are  obtained  allowing  for  the  temperature  dependence  of  viscosity  and  the  dashed  lines  with  a  constant 
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viscosity,  fi  =  fi o-  Clearly  there  is  a  very  large  difference  between  the  two  cases  that  indicates  the  total 
unreliability  of  the  constant-viscosity  model  in  the  later  stages  of  the  process.  Constant  viscosity  results  will 
not  be  considered  further.  Figure  6  shows  the  radial  velocity  in  the  plane  of  symmetry  (solid  lines)  and  the 
maximum  temperature  (occurring  at  the  wall  for  r  =  R,  dashed  lines)  as  functions  of  h(t)/h(0).  The  upper 
pair  of  lines  is  for  V0  =  1  m/s,  and  the  lower  one  for  Vo  =  0.5  m/s.  The  solid  curves  for  the  velocity  exhibit 
a  small  local  maximum  in  the  later  stages  of  the  process.  This  is  due  to  the  formation  of  the  shear-band-like 
structure  near  the  wall  that  has  the  effect  of  causing  a  temporary  slowing  down  of  the  liquid  at  the  midplane. 

Figures  7  to  9  refer  to  the  isothermal  wall  case.  Figure  7  shows  the  radial  velocity  profiles  at  r  =  R 
at  different  times  for  Vo  =  1  m/s.  Qualitatively  the  trend  is  quite  similar  to  that  shown  in  Fig.  2  for  the 
adiabatic  wall  case.  However,  the  shear-band-like  structure  forms  somewhat  later,  as  could  be  anticipated. 
The  corresponding  temperature  distributions,  to  be  compared  with  those  of  Fig.  3,  are  shown  in  Fig.  8. 
Because  of  the  relatively  small  thermal  conductivity  of  the  fluid,  there  is  a  substantial  temperature  increase 
in  this  case  as  well,  although  the  magnitude  is  about  one  half  of  that  found  in  the  adiabatic-wall  case.  The 
velocity  (solid  lines)  and  temperature  histories  comparable  to  those  of  Fig.  6  are  shown  in  Fig.  9.  Upon 
comparison  of  the  two,  it  is  is  seen  that  the  maximum  velocities  are  not  very  strongly  affected  by  the  wall 
boundary  condition,  while  the  temperature  increases  are  reduced  by  about  50%. 

We  now  turn  to  the  more  realistic  situation  in  which  the  motion  of  the  plates  is  opposed  by  the  pressure 
field  developing  in  the  liquid.  We  use  the  same  parameter  values  as  before,  but  now  Vo  is  the  initial  condition 
for  the  velocity  V(t)  that  is  calculated  as  explained  before.  The  major  difference  with  the  previous  case  is 
that  now  the  velocity  of  the  plates  starts  decreasing  appreciably  as  the  pressure  builds  up.  Of  course,  this 
effect  sets  in  later  and  later  as  the  inertia  M  of  the  plates  is  increased,  but  nevertheless  it  is  a  major  effect 
for  any  M.  Velocity  and  temperature  distributions  are  similar  to  those  shown  in  the  previous  figures,  except 
for  the  fact  that  it  is  difficult  to  generate  shear-band  structures.  To  appreciate  the  differences  with  the 
previous  cases  it  will  be  sufficient  to  consider  the  time  histories  of  the  radial  velocity  at  z  =  0  (solid  lines) 
and  maximum  temperature  shown  in  Figs.  10  and  11  for  the  adiabatic  and  isothermal  cases  respectively. 
For  all  cases  Vb  =  1  m/s  and  the  lines  are,  in  descending  order,  for  M  =  2,  1  and  0.2  kg.  A  very  striking 
difference  with  the  previous  situation  is  the  drastic  velocity  reduction.  Naturally,  this  is  accompanied  by 
a  strong  reduction  in  the  maximum  temperature  increase.  The  50%  difference  between  the  adiabatic  and 
isothermal  wall  cases  found  before  is  also  encountered  here. 

It  is  interesting  to  note  that  the  acceleration  becomes  strongly  negative  in  the  later  stages  of  the  flow.  In 
a  situation  in  which  a  single  drop  of  liquid  is  squeezed  between  the  plates,  the  free  surface  along  the  rim  of 
the  drop  would  thus  be  subject  to  an  acceleration  directed  from  the  surrounding  air  into  the  liquid,  which 
is  a  Rayleigh-Taylor  unstable  situation.  Thus  one  would  expect  a  roughening  of  the  drop  periphery  as  is 
indeed  observed  experimentally  (see  e.g.  Field  et  al.  1982).  As  a  matter  of  fact,  the  ignition  events  reported 
by  Field  and  co-workers  all  seem  to  originate  in  this  peripheral  region,  possibly  due  to  air  entrapped  by  the 
surface  deformation  of  the  drop. 

The  final  aspect  that  we  consider  is  that  of  wall  slip.  Our  justification  is  twofold.  In  the  first  place, 
the  usual  no-slip  condition  is  an  empirical  fact  that  has  never  been  tested  at  the  extremely  high  velocities 
encountered  in  some  of  the  cases  that  we  have  simulated.  Secondly,  it  is  known  (see  e.g.  Richardson  1973; 
Davis  and  Miksis  1994)  that  the  no-slip  condition  on  a  rough  wall  may  be  approximately  replaced  by  a  partial 
slip  condition  of  the  form  (3.2.7)  over  a  plane  surface.  In  the  first  case,  the  parameter  A  may  be  expected 
to  be  of  the  order  of  molecular  (or  molecular  clusters)  dimensions,  so  that  A  ~  1-10  nm.  In  the  second  one, 

A  is  of  the  order  of  the  total  volume  of  the  roughnesses  per  unit  area  of  the  surface.  For  a  polished  metal 
surface  typical  values  are  in  the  range  of  0.1  to  1  fim. 

With  slip,  there  is  intense  energy  dissipation  at  the  wall  and  it  would  be  meaningless  to  consider  the  case 
of  an  isothermal  wall.  Hence  we  show  in  Fig.  12  the  wall  temperature  at  r  =  R,  z  =  h(t)  as  a  function 
of  h(t)/h0  for  A  =  100  nm,  10  nm,  and  0.  The  10  nm  case  exhibits  only  a  negligible  difference  with  the 
previous  no-slip  results.  For  A  =  100  nm,  on  the  other  hand,  the  temperature  rise  is  approximately  double 
that  found  in  the  no-slip  case. 
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3.4  Conclusions  on  Viscous  Heating 

In  this  section  we  have  studied  the  viscous  heating  of  a  thin  (~  2  mm)  liquid  layer  squeezed  between  two 
rapidly  approaching  parallel  rigid  plates.  The  main  conclusions  are  as  follows: 

1.  Accounting  for  the  temperature  dependence  of  viscosity  is  essential.  Constant-viscosity  simulations 
give  maximum  temperature  rises  that  can  be  an  order  of  magnitude  larger  than  those  found  allowing 
viscosity  to  decrease  with  increasing  temperature. 

2.  If  the  velocity  of  the  approaching  plates  is  held  constant,  large  temperature  rises  can  be  found  by  the 
time  the  plate  spacing  is  of  the  order  of  10  /im.  This  happens  with  both  adiabatic  and  isothermal  wall 
boundary  conditions,  although  maximum  temperature  rises  are  about  50%  smaller  in  the  latter  case. 
In  addition,  high-velocity  “jets”  resembling  shear  bands  in  solids  form  near  the  boundary  due  to  the 
strong  viscosity  reduction  in  this  region  of  high  shear. 

3.  If  the  reaction  of  the  liquid  on  the  plates  through  the  developing  pressure  field  is  accounted  for,  the 
violence  of  the  phenomenon  is  greatly  reduced  and  typical  maximum  temperature  rises  are  of  the  order 
of  300  0  for  an  adiabatic  wall  and  half  as  much  for  an  isothermal  wall. 

4.  The  wall  temperature  is  found  to  be  very  strongly  sensitive  to  wall  slip  above  a  certain  level.  This 
aspect  of  the  problem  could  be  significant,  but  it  has  not  been  possible  to  pursue  it  to  any  great  depth 
in  view  of  the  fact  that  very  little  seems  to  be  known  about  the  magnitude  of  the  wall  slip  coefficient. 
In  this  connection,  we  might  refer  again  to  the  ignition  of  propellants  near  the  rim  of  the  expanding 
drop  observed  in  experiments.  Indeed,  it  seems  to  be  fairly  well  established  that  the  motion  of  an 
air-liquid-solid  contact  line  is  accompanied  by  local  slip  (see  e.g.  Huh  and  Scriven  1971;  Richardson 
1973;  Dussan  V.  1979). 

These  results  indicate  that  the  process  under  consideration  is  strongly  influenced  by  subtle  aspects  that 
cannot  be  ignored  even  as  a  first  approximation.  In  particular,  the  effect  of  the  large  pressures  developed 
in  the  liquid  has  been  found  to  be  important.  In  the  model  that  we  have  studied,  the  solid  surfaces  were 
assumed  to  be  rigid  and  the  only  effect  of  these  pressures  was  to  slow  down  the  approach  of  the  plates  at  a 
rate  dependent  on  their  inertia.  More  generally,  when  the  inertia  is  very  large,  the  high  local  pressures  will 
lead  to  a  deformation  of  the  approaching  surfaces  which  will  also  tend  to  strongly  limit  the  heating  of  the 
liquid. 

In  all  the  cases  we  have  studied,  substantial  liquid  heating  only  occurs  when  the  film  has  thinned  to  about 
10  ftm  or  less.  At  speeds  of  1  to  0.1  m/s,  the  gap  would  close  completely  in  0.01  to  0.1  ms.  According  to  the 
estimates  of  Field  (1992;  Field  et  al.  1992),  over  time  and  spatial  scales  of  this  magnitudes,  temperatures  of 
the  order  of  700  K  are  needed  for  ignition.  Our  results  imply  that  such  conditions  are  very  hard  to  realize 
by  relying  solely  on  viscous  heating. 
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Fig.  2  Distribution  of  the  radial  velocity  u  across  the  gap  at  tire  edge  of  the  plates  for  a  constant  velocity  of 
approach  Vq  =  1  m/s  and  adiabatic  walls.  The  lines  are,  in  descending  order,  for  h/h0  =  0.00369,  0.00451, 
0.00551,0.00822,0.0123,0.0223,0.0497,0.819. 
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Fig.  3  Temperature  distribution  across  the  gap  at  the  edge  of  the  plates  for  a  constant  velocity  of  approach 
Vo  =  1  m/s  and  an  adiabatic  wall.  The  lines  are,  in  descending  order,  for  li/h0  =  0.00369,  0.00451,  0.00551, 
0.00822,  0.0123,  0.0223,  0.0497,  0.819. 
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ig.  11  Same  as  Fig.  10  for  isothermal  walls. 
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